//////////////////////////////////////////////////////////////////////////
/// Copyright (c) 2011, Sidorenko 'PC' Alexander                       ///
/// mailto: sidorenko /*dot*/ alexander (at) gmail [dot] com           ///
/// This file is a part of the Radiating Gas Dynamics project          ///
//////////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "SodProblem2d.h"
#include "ColorBar.h"

#include "PlotHelper.h"

const gdScalar T0 = gdScalar(2273.0);
const gdScalar T1 = gdScalar(293.0);

const gdScalar p0 = gdScalar(1.01e5);
const gdScalar p1 = p0 * gdScalar(11.0);

const gdScalar r0 = p0 / T0 / R;
const gdScalar r1 = p1 / T1 / R;

SodProblem2d::SodProblem2d(gdScalar gamma, gdScalar cfl, unsigned int numThreads, BOOL saveFull)
: ProblemN<2>(SOD_2D, gamma, cfl, numThreads, saveFull)
, m_xMax(0.0)
{
}

SodProblem2d::~SodProblem2d()
{

}

void SodProblem2d::initGrid(gdScalar xMax, unsigned int numCells,
							 gdScalar tolerance, unsigned int stepsThreshold)
{
	m_dimensions[0] = numCells;
	m_dimensions[1] = numCells;
	m_xMax = xMax;
	gdScalar dx = xMax / gdScalar(numCells);
	// numCells is a dimension. Here's a real number of cells:
	unsigned int fullNumCells = m_dimensions[0] * m_dimensions[1];
	unsigned int numFaces = m_dimensions[0] * (m_dimensions[1] - 1) +
							m_dimensions[1] * (m_dimensions[0] - 1);
	GodunovSolverMT<2>::initModel(fullNumCells, numFaces, dx, tolerance, stepsThreshold);
}

unsigned int SodProblem2d::getNumPlotTypes() const
{
	return 5;
}

std::string SodProblem2d::getPlotTypeName(unsigned int id) const
{
	static const std::string names[] = { "Density", "Pressure", "Temperature", "Velocity", "Momentum" };
	_ASSERT(id < 5);
	return names[id];
}

void SodProblem2d::plotSolution(CDC& dc, const CRect& dstRect, const CColorBar& clr, unsigned int timeStep, unsigned int plotType)
{
	const TimeStep<2> *ts = NULL;
	if (!getData(timeStep, &ts))
	{
		dc.TextOutA(0, 0, "Data is not ready yet!");
		return;
	}
	_ASSERT(NULL != ts);

	const ProblemN<2>::Grid &gr = ts->grid;

	const gdScalar* start = NULL;
	unsigned int offset = sizeof(ProblemN<2>::Cell) / sizeof(gdScalar);
	gdScalar minVal, maxVal;

	switch(plotType)
	{
		case 0: // Density
			start = &gr[0].r;
			minVal = m_minValues.r;
			maxVal = m_maxValues.r;
		break;

		case 1: // Pressure
			start = &gr[0].p;
			minVal = m_minValues.p;
			maxVal = m_maxValues.p;
		break;

		case 2: // Temterature
			start = &gr[0].T;
			minVal = m_minValues.T;
			maxVal = m_maxValues.T;
		break;

		case 3: // Velocity
			start = &gr[0].v[0];
			minVal = m_minValues.v[0];
			maxVal = m_maxValues.v[0];
		break;

		case 4: // Momentum
			start = &gr[0].u[0];
			minVal = m_minValues.u[0];
			maxVal = m_maxValues.u[0];
		break;

		default:
			_ASSERT(0);
		break;
	}

	CRect tmp = dstRect;
	drawPlot(dc, tmp, clr, 0, m_minCellSize * gdScalar(m_grid.size()), minVal, maxVal);

	gdScalar val = *start;
	gdScalar clrFactor = (val - minVal) / (maxVal - minVal);
	int y = tmp.bottom + int(gdScalar(tmp.top - tmp.bottom) * clrFactor);
	dc.SelectObject(GetStockObject(DC_PEN));

	//dc.SetDCPenColor(RGB(0, 0, 0));
	//dc.MoveTo(tmp.left, tmp.top);
	//dc.LineTo(tmp.right, tmp.top);
	//dc.LineTo(tmp.right, tmp.bottom);
	//dc.LineTo(tmp.left, tmp.bottom);
	//dc.LineTo(tmp.left, tmp.top);

	//dc.MoveTo(tmp.left, y);

	gdScalar fW = gdScalar(tmp.Width() - 1);
	gdScalar fH = gdScalar(tmp.Height() - 1);
	gdScalar fD0 = gdScalar(m_dimensions[0] - 1);
	gdScalar fD1 = gdScalar(m_dimensions[1] - 1);

	COLORREF c = 0;

	resizeDrawBuffer(tmp.Size());
	BYTE* dibData = (BYTE*)m_buffer.getData();

	for (unsigned int j = tmp.bottom - 1; j >= tmp.top; --j)
	{
		for (unsigned int i = tmp.left; i < tmp.right; ++i)
		{
			gdScalar factori = gdScalar(i - tmp.left) / fW;
			gdScalar factorj = gdScalar(j - tmp.top) / fH;
			unsigned int ii = unsigned int(factori * fD0);
			unsigned int jj = unsigned int(factorj * fD1);
			val = start[(ii + jj * m_dimensions[0]) * offset];
			clrFactor = (val - minVal) / (maxVal - minVal);
			c = clr.GetColor(clrFactor);
			*dibData++ = GetBValue(c);
			*dibData++ = GetGValue(c);
			*dibData++ = GetRValue(c);
			dibData++;
//			dc.SetPixel(i, j, c);
			//dc.SetDCPenColor(clr.GetColor(clrFactor));
			//dc.LineTo(x, y);
		}
	}
	m_buffer.draw(&dc, CRect(CPoint(0, 0), m_buffer.getSize()), tmp);
	//const int w = 2;
	//for (int i = dstRect.left; i < dstRect.right; i += w)
	//{
	//	gdScalar factor = gdScalar(i) / gdScalar(dstRect.Width() - 1);
	//	int id = int(factor * gdScalar(gridSize - 1));
	//	gdScalar val = start[id * offset];
	//	gdScalar clrFactor = (val - minVal) / (maxVal - minVal);
	//	dc.FillSolidRect(i, dstRect.top, i + w, dstRect.bottom, clr.GetColor(clrFactor));
	//}
}

void SodProblem2d::generateGrid()
{
	gdScalar dx = m_minCellSize,
			 dy = m_minCellSize;
	unsigned int i, j;

	// Build up cells
	for (j = 0; j < m_dimensions[1]; ++j)
	{
		for (i = 0; i < m_dimensions[0]; ++i)
		{
			GodunovSolverMT<2>::Cell &c = m_grid[i + j * m_dimensions[0]];
			c.dv = dx * dy;
		}
	}
	
	// Build up faces
	unsigned int k = 0;

	// Horizontal
	for (j = 0; j < m_dimensions[1]; ++j)
	{
		for (i = 0; i < m_dimensions[0] - 1; ++i)
		{
			_ASSERT(k < m_faces.size());
			GodunovSolverMT<2>::Face &f = m_faces[k++];
			f.da = dy;
			f.n[0] = gdScalar(1.0);
			f.n[1] = gdScalar(0.0);
			f.fc[0] = i + j * m_dimensions[0];
			f.fc[1] = i + 1 + j * m_dimensions[0];
		}
	}

	// Vertical
	for (j = 0; j < m_dimensions[1] - 1; ++j)
	{
		for (i = 0; i < m_dimensions[0]; ++i)
		{
			_ASSERT(k < m_faces.size());
			GodunovSolverMT<2>::Face &f = m_faces[k++];
			f.da = dx;
			f.n[0] = gdScalar(0.0);
			f.n[1] = gdScalar(1.0);
			f.fc[0] = i + j * m_dimensions[0];
			f.fc[1] = i + (j + 1) * m_dimensions[0];
		}
	}

	_ASSERT(k == m_faces.size());
}

void SodProblem2d::generateInitialDistributions()
{
	unsigned int i, j;


	for (i = 0; i < m_dimensions[0]; ++i)
	{
		for (j = 0; j < m_dimensions[1]; ++j)
		{
			GodunovSolverMT<2>::Cell &c = m_grid[i + j * m_dimensions[0]];
			c.p = 0.1f;
			c.r = 0.125f;
			c.v.set(gdScalar(0.0));
			c.u.set(gdScalar(0.0));
		}
	}

	for (i = 0; i < m_dimensions[0]; ++i)
	{
		for (j = 0; j < i; ++j)
		{
			GodunovSolverMT<2>::Cell &c = m_grid[i + j * m_dimensions[0]];
			c.p = 1.0f;
			c.r = 1.0f;
			c.v.set(gdScalar(0.0));
			c.u.set(gdScalar(0.0));
		}
	}
}

void SodProblem2d::applyBoundaryConditions()
{
	unsigned int i, j;

	for (i = 0; i < m_dimensions[0]; ++i)
	{
		GodunovSolverMT<2>::Cell &c1 = m_grid[i + 0 * m_dimensions[0]];
		GodunovSolverMT<2>::Cell &c2 = m_grid[i + 1 * m_dimensions[0]];
		//c1.r = 0;
		//c1.p = 0;
		//c1.u[0] = 0;
		//c1.u[1] = 0;
		c1.r = c2.r;
		c1.p = c2.p;
		c1.u[0] = c2.u[0];
		c1.u[1] = -c2.u[1];

		GodunovSolverMT<2>::Cell &c3 = m_grid[i + (m_dimensions[1] - 2) * m_dimensions[0]];
		GodunovSolverMT<2>::Cell &c4 = m_grid[i + (m_dimensions[1] - 1) * m_dimensions[0]];
		c4.r = c3.r;
		c4.p = c3.p;
		c4.u[0] = c3.u[0];
		c4.u[1] = -c3.u[1];
	}

	for (j = 0; j < m_dimensions[1]; ++j)
	{
		GodunovSolverMT<2>::Cell &c1 = m_grid[0 + j * m_dimensions[0]];
		GodunovSolverMT<2>::Cell &c2 = m_grid[1 + j * m_dimensions[0]];
		c1.r = c2.r;
		c1.p = c2.p;
		c1.u[0] = -c2.u[0];
		c1.u[1] = c2.u[1];

		GodunovSolverMT<2>::Cell &c3 = m_grid[m_dimensions[0] - 2 + j * m_dimensions[0]];
		GodunovSolverMT<2>::Cell &c4 = m_grid[m_dimensions[0] - 1 + j * m_dimensions[0]];
		c4.r = c3.r;
		c4.p = c3.p;
		c4.u[0] = -c3.u[0];
		c4.u[1] = c3.u[1];
	}
}

void SodProblem2d::solveRadiation(gdScalar timeDelta)
{
	// Here you should solve radiation problem; All
	// the variables in the m_grid, such as Density r,
	// Energy e, Temperature T are up-to-date and valid
	// I GUARANTEE ;)
}